###########################################################################################
###########################################################################################
######################-~-Interest Representation in Belgium-~-#############################
######################     Politics of the Low Countries      #############################
###########################################################################################
#### Evelien Willems ######################################################################
#### Frederik Heylen ######################################################################
#### Jan Beyers      ######################################################################
###########################################################################################
###########################################################################################

#_____________________________________Load packages_______________________________________#
library("ggplot2")
library("psych")
library("dyplr")
library("car")

#__________________________Load data - SET WORKING DIRECTORY______________________________#
  #only the survey data + website coding
  cigBE <- read.csv("H:/3_5_population and CIG-survey project/2_publications/9_Low Countries 2019/replication files/cigBE_PLC.csv")
  
  #with access to advisory councils
  cigBE_AC <- read.csv("H:/3_5_population and CIG-survey project/2_publications/9_Low Countries 2019/replication files/cigBE_AC_PLC.csv")
  
  #with media attention
  cigBE_NM <- read.csv("H:/3_5_population and CIG-survey project/2_publications/9_Low Countries 2019/replication files/cigBE_NM_PLC.csv")

#__________________________Descriptive results_Main text__________________________________#
  
#Table 2 - Demography by group type across government levels
  #distribution across the different levels of government
  prop.table(table(cigBE$q73r))
  prop.table(table(cigBE$q73r, cigBE$type),1)
  prop.table(table(cigBE$type))
  #labels government levels (q73r): 1=national; 2=Flemish; 3=Francophone
  #labels type: 1=business; 2=professional; 3=labour; 4=identity; 5=cause; 6=leisure; 7=institutions
  
#Table 3 - Basic organisational features by organisation type (survey responses only)
  ## Number of groups
  table(cigBE$type & cigBE$respid)
  as.data.frame(table(cigBE$type, cigBE$respid >= 1))
  
  ## Mean Year of founding (q02r)
  mean(na.omit(cigBE[which(cigBE$type==1), which(colnames(cigBE)=="q02r" )]))
  mean(na.omit(cigBE[which(cigBE$type==2), which(colnames(cigBE)=="q02r" )]))
  mean(na.omit(cigBE[which(cigBE$type==3), which(colnames(cigBE)=="q02r" )]))
  mean(na.omit(cigBE[which(cigBE$type==4), which(colnames(cigBE)=="q02r" )]))
  mean(na.omit(cigBE[which(cigBE$type==5), which(colnames(cigBE)=="q02r" )]))
  mean(na.omit(cigBE[which(cigBE$type==6), which(colnames(cigBE)=="q02r" )]))
  mean(na.omit(cigBE[which(cigBE$type==7), which(colnames(cigBE)=="q02r" )]))
  
  ## Mean FTE (q21_01r)
  cigBE$q21_01r <- as.numeric(cigBE$q21_01r)
  mean(na.omit(cigBE[which(cigBE$type==1), which(colnames(cigBE)=="q21_01r" )]))
  mean(na.omit(cigBE[which(cigBE$type==2), which(colnames(cigBE)=="q21_01r" )]))
  mean(na.omit(cigBE[which(cigBE$type==3), which(colnames(cigBE)=="q21_01r" )]))
  mean(na.omit(cigBE[which(cigBE$type==4), which(colnames(cigBE)=="q21_01r" )]))
  mean(na.omit(cigBE[which(cigBE$type==5), which(colnames(cigBE)=="q21_01r" )]))
  mean(na.omit(cigBE[which(cigBE$type==6), which(colnames(cigBE)=="q21_01r" )]))
  mean(na.omit(cigBE[which(cigBE$type==7), which(colnames(cigBE)=="q21_01r" )]))
  
  ## Median Budget (q08r)
  median(na.omit(cigBE[which(cigBE$type==1), which(colnames(cigBE)=="q08r" )]))
  median(na.omit(cigBE[which(cigBE$type==2), which(colnames(cigBE)=="q08r" )]))
  median(na.omit(cigBE[which(cigBE$type==3), which(colnames(cigBE)=="q08r" )]))
  median(na.omit(cigBE[which(cigBE$type==4), which(colnames(cigBE)=="q08r" )]))
  median(na.omit(cigBE[which(cigBE$type==5), which(colnames(cigBE)=="q08r" )]))
  median(na.omit(cigBE[which(cigBE$type==6), which(colnames(cigBE)=="q08r" )]))
  median(na.omit(cigBE[which(cigBE$type==7), which(colnames(cigBE)=="q08r" )]))
  #labels budget categories (q08r): 1= >10.000; 2=10.000-50.000; 3=50.000-100.000; 4=100.000-500.000; 5=500.000-1.000.000; 6=1.000.000-5.000.000; 7=5.000.000-10.000.000; 8= <10.000.000
  
  
  
#Figure 1 - Founding dates of Belgian interest organisations
  ggplot(cigBE_AC, aes(YoF , na.rm = TRUE)) + geom_histogram(breaks = seq(1830, 2020, by = 5), col="black", fill="grey", alpha = 1, binwidth = 5,  na.rm =TRUE )+
    labs(title="") + 
    labs(x="Year", y="Number of founded organisations") +
    scale_x_continuous(breaks=seq(1830, 2020, 10)) +
    theme(axis.title.y = element_text(color="black", size=14, face="bold"),
          axis.title.x = element_text(color="black", size=14, face="bold"))+
    theme_classic()
  
#Figure 2 - Inside and outside advocacy tactics (bar graph created in excel, category 4 and 5 taken together)
  #direct contacts
  table(cigBE$directcontacts)
  #stakeholder meeting
  table(cigBE$q34_10r)
  #website: position papers
  table(cigBE$q34_08r)
  #protest
  table(cigBE$q34_07r)
  #member action
  table(cigBE$q34_06r)
  #contact journalists
  table(cigBE$q34_05r)
  #media appearances
  table(cigBE$q34_03r)
  #research reports
  table(cigBE$q34_02r)
  #press conferences
  table(cigBE$q34_01r)
  #together
  tactics <- rbind(table(cigBE$directcontacts), table(cigBE$q34_10r), table(cigBE$q34_08r), table(cigBE$q34_07r),table(cigBE$q34_06r),
                   table(cigBE$q34_05r), table(cigBE$q34_03r), table(cigBE$q34_02r), table(cigBE$q34_01r))
  tactics
  #labels frequency of tactics used: 1=never; 2=at least once a year; 3=once every 3 months; 4=once every month; 5=once per week
  
#Figure 3 - Seeking access to different governmental branches by governmental level (bar graph created in excel, category 3 and 4 taken together) 
  #labels frequency of contact-seeking: 1=never; 2=at least once a year; 3=once every 3 months; 4=once every month; 5=once per week
  #executive national
  a <- as.data.frame(prop.table(table(cigBE$q33_01rr))) # ministers themselves
  b <- as.data.frame(prop.table(table(cigBE$q33_04rr))) # office of the prime minister 
  c <- as.data.frame(prop.table(table(cigBE$q33_07rr))) # European affairs ministers' offices
  d <- as.data.frame(prop.table(table(cigBE$q33_08rr))) # other ministers' offices
  
  A <- (a[1,2] + b[1,2] + c[1,2] + d[1,2])/4
  B <- (a[2,2] + b[2,2] + c[2,2] + d[2,2])/4
  C <- (a[3,2] + b[3,2] + c[3,2] + d[3,2])/4
  D <- (a[4,2] + b[4,2] + c[4,2] + d[4,2])/4
  E <- (a[5,2] + b[5,2] + c[5,2] + d[5,2])/4
  
  #mean contact to the executive national
  fedexec <-  c(A,B,C,D,E)
  fedexec
  
  #bureaucracy national
  a <- as.data.frame(prop.table(table(cigBE$q33_05rr))) # civil servants
  b <- as.data.frame(prop.table(table(cigBE$q33_06rr))) # civil servants
  
  A <- (a[1,2] + b[1,2])/2
  B <- (a[2,2] + b[2,2])/2
  C <- (a[3,2] + b[3,2])/2
  D <- (a[4,2] + b[4,2])/2
  E <- (a[5,2] + b[5,2])/2
  
  fedbur <-  c(A,B,C,D,E)
  fedbur
  
  #parliament national
  a <- as.data.frame(prop.table(table(cigBE$q33_02rr))) # MPs governing parties
  b <- as.data.frame(prop.table(table(cigBE$q33_03rr))) # MPs opposition parties
  
  A <- (a[1,2] + b[1,2])/2
  B <- (a[2,2] + b[2,2])/2
  C <- (a[3,2] + b[3,2])/2
  D <- (a[4,2] + b[4,2])/2
  E <- (a[5,2] + b[5,2])/2
  
  fedpar <-  c(A,B,C,D,E)
  fedpar
  
  #executive francophone (Walloon region + Francophone community)
  a <- as.data.frame(prop.table(table(cigBE$q100_01rr))) # ministers themselves
  b <- as.data.frame(prop.table(table(cigBE$q100_02rr))) # office of the prime minister
  c <- as.data.frame(prop.table(table(cigBE$q100_03rr))) # other ministers' offices
  
  aa <- as.data.frame(prop.table(table(cigBE$q101_01rr))) # ministers themselves
  bb <- as.data.frame(prop.table(table(cigBE$q101_02rr))) # office of the prime minister
  cc <- as.data.frame(prop.table(table(cigBE$q101_03rr))) # other ministers' offices
  
  A <- (a[1,2] + b[1,2] + c[1,2] + aa[1,2] + bb[1,2] + cc[1,2])/6
  B <- (a[2,2] + b[2,2] + c[2,2] + aa[2,2] + bb[2,2] + cc[2,2])/6
  C <- (a[3,2] + b[3,2] + c[3,2] + aa[3,2] + bb[3,2] + cc[3,2])/6
  D <- (a[4,2] + b[4,2] + c[4,2] + aa[4,2] + bb[4,2] + cc[4,2])/6
  E <- (a[5,2] + b[5,2] + c[5,2] + aa[5,2] + bb[5,2] + cc[5,2])/6
  
  #mean contact to the executive
  walloonexec <-  c(A,B,C,D,E)
  walloonexec
  
  #bureaucracy francophone (Walloon region + Francophone community)
  a <- as.data.frame(prop.table(table(cigBE$q100_04rr))) # civil servants
  b <- as.data.frame(prop.table(table(cigBE$q101_04rr))) # civil servants
  
  A <- (a[1,2] + b[1,2])/2
  B <- (a[2,2] + b[2,2])/2
  C <- (a[3,2] + b[3,2])/2
  D <- (a[4,2] + b[4,2])/2
  E <- (a[5,2] + b[5,2])/2
  
  walloonbur <-  c(A,B,C,D,E)
  walloonbur
  
  #parliament francophone (Walloon region + Francophone community)
  a <- as.data.frame(prop.table(table(cigBE$q100_05rr))) # MPs governing parties
  b <- as.data.frame(prop.table(table(cigBE$q100_06rr))) # MPs opposition parties
  c <- as.data.frame(prop.table(table(cigBE$q101_05rr))) # MPs governing parties
  D <- as.data.frame(prop.table(table(cigBE$q101_06rr))) # MPs opposition parties
  
  A <- (a[1,2] + b[1,2] + c[1,2] + d[1,2])/4
  B <- (a[2,2] + b[2,2] + c[2,2] + d[2,2])/4
  C <- (a[3,2] + b[3,2] + c[3,2] + d[3,2])/4
  D <- (a[4,2] + b[4,2] + c[4,2] + d[4,2])/4
  E <- (a[5,2] + b[5,2] + c[5,2] + d[5,2])/4
  
  walloonpar <-  c(A,B,C,D,E)
  walloonpar
  
  #executive flemish
  a <- as.data.frame(prop.table(table(cigBE$q103_01rr))) # ministers themselves
  b <- as.data.frame(prop.table(table(cigBE$q103_02rr))) # office of the prime minister
  c <- as.data.frame(prop.table(table(cigBE$q103_03rr))) # other ministers' offices 
  
  
  A <- (a[1,2] + b[1,2] + c[1,2])/3
  B <- (a[2,2] + b[2,2] + c[2,2])/3
  C <- (a[3,2] + b[3,2] + c[3,2])/3
  D <- (a[4,2] + b[4,2] + c[4,2])/3
  E <- (a[5,2] + b[5,2] + c[5,2])/3
  
  #mean contact to the executive
  flexec <-  c(A,B,C,D,E)
  flexec
  
  #bureaucracy flemish
  flbur <- prop.table(table(cigBE$q103_04rr)) # civil servants
  flbur
  
  #parliament flemish
  a <- as.data.frame(prop.table(table(cigBE$q103_05rr))) # MPs governing parties
  b <- as.data.frame(prop.table(table(cigBE$q103_06rr))) # MPs opposition parties
  
  A <- (a[1,2] + b[1,2])/2
  B <- (a[2,2] + b[2,2])/2
  C <- (a[3,2] + b[3,2])/2
  D <- (a[4,2] + b[4,2])/2
  E <- (a[5,2] + b[5,2])/2
  
  flpar <-  c(A,B,C,D,E)
  flpar
  
  contacts <- rbind(fedexec, fedbur, fedpar, walloonexec, walloonbur, walloonpar, flexec, flbur, flpar)
  contacts 
  
#Figure 4 - Comparing access of different group types to advisory councils (bar graph created in excel)
  #Use data-file 'cigBE_AC'
  cigBE_AC$seats <- cigBE_AC$X..E + cigBE_AC$X..A + cigBE_AC$X..O  #summing all types of seats per organisation: effective, alternate, observer                     
  
  cigBE_AC <- within (cigBE_AC,{
    seats1 <- Recode(seats, '0 = NA'
    )}) #Exclusion of interest groups without access to advisory councils
  
  describeBy(cigBE_AC$seats1, cigBE_AC$type3, quant=c(0.25, 0.75))
  table(cigBE_AC$type3 & cigBE_AC$seats)
  table(cigBE_AC$type3)
  as.data.frame(table(cigBE_AC$type3, cigBE_AC$seats > 0))
  
#Table 4 - List of top 20 advisers
  #Script not included to protect anonymity of survey respondents
  
#Figure 5 - Comparing media attention for different group types 
  #Use data-file 'cigBE_NM'
  describeBy(cigBE_NM$accessmedia, cigBE_NM$type3, quant=c(0.25, 0.75))
  table(cigBE_NM$type3 & cigBE_NM$accessmedia)
  table(cigBE_NM$type3)
  as.data.frame(table(cigBE_NM$type3, cigBE_NM$accessmedia > 0))
  
#Table 5 - The distribution of media attention by group type, at least one media hit   
  describeBy(cigBE_NM$accessmedia, cigBE_NM$type3, quant=c(0.25, 0.75))
  
#________________________________________Descriptive results_Online Appendix__________________________________________#
#Figure 2A - Founding dates by group type
  YOF1 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==1),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Business associations", 
                xlab = "",  
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  YOF2 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==2),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Professional associations", 
                xlab = "",
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  YOF3 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==3),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Labour unions", 
                xlab = "", 
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  YOF4 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==4),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Identity groups", 
                xlab = "", 
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  YOF5 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==5),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Cause groups", 
                xlab = "", 
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  YOF6 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==6),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Leisure associations", 
                xlab = "",
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  YOF7 <- qplot(na.omit(cigBE[which(cigBE$YoF > 1830 & cigBE$type ==7),which(colnames(cigBE)=="q02r" )]),
                geom="histogram",
                binwidth = 5, 
                main = "Founding dates Institutions associations", 
                xlab = "",
                ylab= "",
                fill=I("grey"), 
                col=I("black"),
                alpha=I(.5)) +
    theme_light() +
    scale_x_continuous(breaks=seq(1830, 2020, 20))
  
  library(gridExtra)
  grid.arrange(YOF1, YOF2, YOF3, YOF4, YOF5, YOF6, YOF7)
  
  
  